function y = func_Sm(zeta, m, CONSTS)

    debug_output = false;
    
    k0 = CONSTS.k0;
    a  = CONSTS.a;
    d  = CONSTS.d;

    y = zeros(size(zeta));
    
    for i = 1:size(zeta, 1)
        zeta_i = abs(zeta(i));

        period = pi/(k0*a);
        spans_per_period = 8;
        order = 1; %3
        %q_high = 0.5*log(10)*order/(k0*zeta_i + 0.001*d);
        if (size(zeta, 1) >= 2)
            adding = 0.9*k0*zeta(2);
        else
            adding = 0.05*d;
        end
        q_high = 0.5*log(10)*order/(k0*(zeta_i) + adding);

        span_size = period/spans_per_period;
        
        %integ_step = 0.05;
        %q_head = (0.0+fallback : integ_step : 100)';
        %q_tail = (100 : integ_step : q_high)';
        
        %integ_head = trapz(q_head, imag(func_ns(q_head, zeta_i, m, CONSTS)));

        integ = quad_gauss(@(q)(func_ns(q, zeta_i, m, CONSTS)), ...
            [0.0, q_high], span_size, 4);
        %integ_tail = trapz(q_tail, imag(func_ns(q_tail, zeta_i, m, CONSTS)));

        if (debug_output)
            fprintf('i=%d    %e + %ei\n', i, ...
                real(integ), imag(integ));
        end
        
        %integ_full = integ_head + integ_tail;

        y(i) = integ;
        
    end
    
end